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ABSTRACT 

In  principle,  given  the  potential  energy  function,  the  values  of  thermo¬ 
dynamic  variables  can  be  computed  from  statistical  mechanics  for  a  system  of 
molecules.  hi  practice  for  the  liquid  state,  however,  two  barriers  must  be  over- 
-^esieX  This  paper  treats  the  first  problem,-  how  to  quantum  correct  the  classical 
mechanical  thermodynamic  values  available  from  molecular  dynamics,  Monte 
Carlo,  perturbation,  or  integral  methods  in  order  to  compare  with  experimental 
quantum  reality.  A  subsequent  paper  will  focus  on  the  second  difficulty,  the 
effective  computation  of  free  energy  and  entropy.  A  simple  technique,  derived 
from  spectral  analysis  of  the  atomic  velocity  time  histories,  is  presented  here 
Jor  the  quantum  correction  of  classical  thermodynamic  values.  This  technique 
ris  based,  on  the  approximation  that  potential  anharmonirities  mainly  affect  the 
lower  frequencies  in  the  velocity  spectrum  where  the  system  behaves  essentially 
'Cigssically,  while  the  higher  spectral  frequencies,  where  the  deviation  from  das- 
step  mechanics  is  most  pronounced,  involve  sufficiently  harmonic  atomic 
mdtioos  that  harmonic  quantum  corrections  apply.  The  approach  is  demon¬ 
strated  by  computation  of  the  energy  and  constant  volume  heat  capacity  for 
water  from  classical  molecular  dynamics  followed  by  quantum  correction.  The 
potential  used  to  describe  the  interactions  of  the  system  of  water  molecules 
indudes  internal  vibrational  degrees  of  freedom  and  thus  strong  quantum 
effects.  Comperison  of  the  quantum  corrected  theoretical  values  with  experi¬ 
mental  measurements  shows  good  agreemenu/fae  quantum  eorrections  to 
classical  thermodynamics  (which  are  also  deri-wto  free  energy  and  entropy) 
are  shown  to  be  important  not  only  for  internal  vibrational  motion,  but  also  for 
intermoiecular  hindered  rotational  and  translational  motions  in  liquid  water. 
They  are  presumably  also  important  for  other  strongly  associated  molecules, 
induding  biomolecules,  and  thus  should  be  induded  when  comparing  calculated 
and  measured  thermodynamic  quantities.  The  approach  illustrated  here  allows 
the  calculation  of  thermodynamic  quantum  corrections  for  liquids,  solutions, 
and  large  molecules  such  as  polymers  (induding  proteins  and  nucleic  adds) 
with  full  indusion  of  both  intra-  and  intermoiecular  degrees  of  freedom. 
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I.  INTRODUCTION 

In  principle  from  the  potential  energy  as  a  function  of  nudear  positions  one  can  compute 
from  statistical  mechanics  the  values  of  the  thermodynamic  variables.  In  practice  this  has  been 
a  difficult  task  for  liquids  and  larger  molecules  such  as  proteins  and  nucleic  adds.  Two  substan¬ 
tial  barriers  need  to  be  overcome.  The  first,  which  is  the  subject  of  this  paper,  is  how  to  com¬ 
pute  quantum  thermodynamic  reality  when  only  classical  mechanics  is  practically  available  as  a 
computational  tool.  That  quantum  mechanics  is  essential  in  treating  intramolecular  vibrations  is 
universally  acknowledged,  but  it  has  sometimes  been  less  well  appreciated  that  intermoiecular 
motions  in  strongly  aesodared  liquids  like  water  also  show  important  quantum  effects.  Quantum 
corrections  should  thus  be  considered  for  strongly  interacting  molecules  in  general,  even  for 
moiecuies  approximated  as  rigid  bodies,  and  for  bkxnotecuies.  The  second  barrier,  which  is  the 
subject  of  a  paper  to  follow,  is  how  to  practically  compute  the  useful,  but  intrihsicaHy  difficult, 
free  energy  and  entropy. 

The  present  paper  illustrates  a  simple  molecular  dynamics  technique  for  quantum  correct¬ 
ing  datsicai  thermodynamic  quantities,  for  example  those  derived  from  molecular  dynamics, 
Monte  Carlo,  perturbation,  or  integral  methods.  This  approach  makes  use  of  the  vdodty  spec¬ 
trum  (often  called  the  velocity  autocorrelation  spectrum),  which  is  related  to  infrared,  Raman, 
and  inelastic  neutron  spectre.  For  harmonic  systems  die  vdodty  spectrum  is  directly  linked  to 
both  classical  and  quantum  mechanical  thermodynamic  parameters,  as  it  then  represents  the 
density  of  normal  mode  harmonic  oscillators  as  a  function  o(  frequency.  Two  suppositions  are 
used  to  justify  a  harmonic  approach  to  estimating  the  thermodynamic  quantum  corrections:  /) 
that  anharmonidttes  mainly  affect  the  low  frequency  motions  which  are  nearly  classical,  and  U) 
that  high  frequency  motions,  where  quantum  effects  are  more  important,  are  nearly  harmonic. 
With  these  assumptions  the  quantum  correction*  frx  a  thermodynamic  variable  can  be 
evaluated  simply  from  the  integral  over  frequen  universal  weighting  function  for  that 

variable  times  the  vdodty  spectrum  computed  fron  spectra  of  atomic  vdodty  time  his¬ 
tories.  The  weighting  ftmcrions  approach  «ro  in  th*.  -  jquency  region  where  anharmonici- 
ties  would  otherwise  cause  problems.  Such  a  quantum  correction  approach  is  not  limited,  like 
most  other  approaches,  to  nearly  classical  systems,  but  can  equally  be  used  to  treat  molecular 
systems  with  internal  vibrational  degrees  of  freedom  where  quantum  effects  are  very  strong,  for 
example  molecular  liquids,  solutions,  solids,  and  polymers,  including  proteins  and  nuddc  adds, 
with  full  indusion  of  internal  degrees  of  freedom. 

Section  Q  describes  the  classical  calculation  of  energy,  heat  capadty,  free  energy,  and 
entropy  from  molecular  dynamics,  followed  in  Section  III  with  the  theory  of  our  quantum 
correction  technique.  Section  IV  describes  the  calculation  and  quantum  correction  of  the 
energy  and  heat  capadty  of  liquid  water.  While  quite  good  agreement  is  achieved  with  experi¬ 
ment,  we  emphasize  that  our  main  purpose  is  to  illustrate  the  techniques  and  not  to  make  the 
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most  accurate  possible  thermodynamic  calculations.  We  point  out  that  the  choice  of  boundary 
treatment  can  significantly  affect  the  numerical  results.  Section  V  discusses  these  results  and 
their  meaning. 

II.  CLASSICAL  THERMODYNAMICS  FROM  MOLECULAR  DYNAMICS 

The  standard  equations1  linking  the  canonical  partition  function  Q  and  the  various  ther¬ 
modynamic  variables  are 

(2.1) 

C.-|f  (2.2) 

A--kgT\nQ  (2.3) 

5  -  kgT~~&  +  kginQ,  (2.4) 

in  which  £  is  the  energy,  C,  the  constant  volume  heat  capacity,  a  the  Helmholtz  free  energy, 
5  the  entropy,  kg  Boitzman’s  constant,  and  T  the  temperature. 

A.  Energy 

The  energy  at  a  given  temperature  T  may  be  computed  from  molecular  dynamics  in 
several  ways.  /)  Choose  initial  conditions  for  a  set  of  different  constant  energy  microcanonical 
molecular  dynamics  runs  to  approximate  a  canonical  ensemble  at  T,  for  example  by  a  sequence 
of  kinetic  energy  randomizations  from  a  Boltzmann  distribution.  The  classical  energy  of  the 
system 

£-<£k(pA0+  K(r*)>,  (2.5) 

is  then  derived  as  an  averaga,  symbolized  by  <>,  over  several  molecular  dynamics  runs  from 
the  ensemble  at  temperature  T  in  which  Eg  is  the  kinetic  energy  and  V  the  potential  energy, 

letting  die  positions  and  momenta  of  the  N  atoms  be  represented  by  rv  m  tt . r*  and 

piV  ■  p], . . .  ,p*.  respectively.  If)  Compute  the  temperature  for  several  different  runs  at 
different  constant  values  of  the  total  energy  by  averaging  the  instantaneous  temperature  defined 
in  terms  of  the  kinetic  energy  by 


T(t)>  -  (3Nkg)-'  £  m,<jv,(r))2>, 


where  kg  is  Boitzman’s  constant,  v,  is  a  Cartesian  component  of  the  velocity  of  one  of  the  N 
moms,  mj  is  the  mass  of  that  atom,  and  <>  here  indicates  a  time  average.  Fit  an  energy 
versus  temperature  curve  to  the  results  for  several  such  microcanonical  molecular  dynamics 
runs.  IB)  Adjust  the  kinetic  energies  during  each  molecular  dynamics  run  in  order  to  represent 
the  system  in  a  heat  bath  at  temperature  T  as  demonstrated  by  Andersen.2  In  this  paper  we  use 
both  approaches  I)  and  B). 

B.  Heat  capacity 

By  performing  microcanonical  molecular  dynamics  rum  at  several  different  energies  and 
computing  the  average  temperature  for  each  energy,  in  other  words  method  //)  above,  the  heat 
capacity  at  constant  volume  C*  can  be  derived  through  numerical  differentiation  of  energy  £ 
with  respect  to  the  temperature  T. 

In  addition,  the  heat  capacity  may  be  calculated  in  principle  from  the  kinetic  energy 
fluctuation  for  a  microcanonical  ensemble.  With  the  velocity  of  the  center  of  mass  set  to  zero, 
the  heat  capacity  is  given  by3*4 


C-R/  -f-tf 


<P>-<T> * 
<T>1 


.*■  "...ini **»»>..■**'  «*  *_#" •* 


a 


in  which  /?  is  kB  times  Avogadro’s  number,  the  number  of  atoms  is  S,  and  T  is  defined  as  in 
Eq.  (2.6)  above.  Statistical  accuracy  becomes  very  important  as  the  denominator  becomes 
small,  which  possibly  explains  why  we  did  not  succeed  in  calculating  accurate  values  using  this 
approach.  . 

C.  Free  energy  and  entropy 

The  free  energy  may  be  computed  from  molecular  dynamics  by  a  technique  due  to  Kirk- 
wood1- 5  which  has  been  applied  in  a  parallel  manner  to  Monte  Carlo  calculations,6  as  demon¬ 
strated  by  Mezei,  Swaminathan  and  Beveridge7-*  in  a  classical  Monte  Carlo  calculation  of  the 
free  energy  for  rigid  molecule  liquid  water. 

The  classical  canonical  ensemble  partition  function  0(f)  is  defined  as1 

<?(f)  -  (Nlh*")-'  J /  drW  expH)//(rVv,f)],  (2.8) 

in  which  H( r‘v,p‘v,f)  is  the  classical  Hamiltonian  of  the  system,  the  Kirkwood1-5*7  f  is  a 
parameter  upon  which  the  Hamiltonian  depends,  and  fi  m  (kaD~l.  Eq.  (2.3)  now  gives 

vt(f)--r'lnO(f).  (2.9) 

Differentiating  Eq.  (19)  with  respect  to  f  gives 


(2.11) 


idifi  _  -0-1  USfilfi  tt» 

which  allows  us  to  write 

A  (f2)  -  A  (f ,)  -  -r‘  ,  (2.11) 

in  which  fj  is  the  value  of  the  Kirkwood  parameter  which  gives  the  real  Hamiltonian  and  ft  is  a 
value  which  distorts  the  Hamiltonian  to  give  a  reference  system  (for  example  an  ideal  gas,  a 
hard  sphere  liquid,  or  a  harmonic  solid)  for  which  we  can  more  easily  compute  the  free 
energy.9  thing  Eq.  (18),  we  have 

mssi _ i  mi 

«f  <ReF  9f 

-fiWi!3*)-'} fdr*dp*  expH3/f(rv,p*,f)l 

- cifl -  ai2) 

which  by  the  ensemble  postulate  of  Gibbs 


P  94 


(113) 


where  the  derivative  of  the  Hamiltonian  with  respect  to  f  is  averaged  over  coordinates  and 
momenta  from  an  ensemble  with  the  Hamiltonian  containing  the  parameter  f.  Substituting  Eq. 
(113)  into  Eq.  (Ill)  gives 


A  (f2)  -  A (#,)  -  jfrff 


(2.14) 


To  evaluste  Eq.  (114)  by  molecular  dynamics,  atomic  trajectories  are  computed  for  the  Hamil¬ 
tonian  /f(rv,pv,f),  &H2SLJL  •(?  is  averaged  over  an  ensemble  of  these  trajectories  at  tem¬ 
perature  r,  and  the  result  is  then  integrated  between  f  j  and  f2. 

In  this  way,  the  classical  free  energy  change  between  the  system  with  our  real  Hamil¬ 
tonian  /f(rv.pv,f2)  and  a  reference  system  with  Hamiltonian  /f(rv.p  v.f ,)  can  be  computed. 


We  chose  the  reference  system  to  be  one  for  which  we  can  compute  the  classical  free  energy 
more  tractably. 

The  entropy  S  of  the  system  may  then  be  calculated  using 
S-(E-A)/T.  (2.16) 

We  will  illustrate  in  a  subsequent  paper  the  actual  molecular  dynamics  calculation  and 
quantum  correction  of  the  free  energy  and  entropy  of  liquid  water  using  approaches  based  on 
the  Kirkwood  technique. 

III.  QUANTUM  CORRECTIONS  FROM  CLASSICAL  MOLECULAR  DYNAMICS 

Outside  of  the  trivial  correction  for  vibrational  zero  point  energy  which  may  be  calculated 
from  spectroscopic  data10  and  which  is  generally  introduced  as  a  constant  in  the  potential  func¬ 
tion,  die  vast  mqjority  of  work  in  quantum  corrections  to  classical  thermodynamic  computations 
sums  from  a  method  first  introduced  by  Wigner11  and  Kirkwood.12- 12  In  this  approach  the  free 
energy  is  expanded  in  powers  of  A2,  and  the  first  term  in  the  quantum  correction  to  be  added  to 
the  rlarairal  value  of  the  free  energy  is  shown  to  be  proportional  to  the  classically  averaged  sum 
of  the  squares  of  the  forces  exerted  on  the  particles  in  the  system.  The  Wigner-Kirkwood  tech¬ 
nique  has  bean  motfified,  extended  and  tested  by  many  workers.14*22  Others22*32  have  examined 
various  methods  to  handle  nondifferentiable  potential  functions  which  apply,  for  example,  to 
hard  spheres  or  square  wells.  Barker  and  Henderson  have  written  a  comprehensive  review  of 
fiquMs  which  includes  an  extensive  section  on  quantum  corrections.6 

Another  quantum  correction  method  by  Doll  and  Myers33  is  based  on  the  path  integral 
approach  of  Feynman  and  Hlbbs.34  It  involves  the  calculation  of  an  effective  potential  Vtff  in 
the  first  stags  of  a  Monte  Carlo  technique.  In  the  second  stage,  y./f  is  used  to  calculate  the 
ratio  between  the  quantum  mechanical  and  classical  partition  functions.  Stillinger35  discusses 
the  sealer  calculation  of  effective  potentials  for  pairwise  potentials. 

h  addition  to  tha  quantum  corrections  considered  here  there  are  the  effects  of  the  sym¬ 
metry  restrictions  on  quantum  states  imposed  by  Fermi-Dirac  and  Boee- Einstein  statistics.  In 
the  temperature  range  of  interest  here  these  effects  are  negligible.12- 13  36 

A  disadvantage  of  all  the  previously  cited  techniques,  except  the  vibrational  zero  point 
energy  correction.  Is  diet  they  are  ordinarily  restricted  to  systems  with  small  quantum  effects. 
The  method  we  present  b  this  paper  may  be  applied  when  quantum  corrections  are  large,  for 
example  to  intramolecular  vibrations. 

Owidd  and  Scheraga37  discuss  the  quantum  corrections  for  liquid  water.  Using  approxi¬ 
mations  to  the  effects  of  Ubrational  and  vibrational  frequencies,  they  calculate  the  quantum 
mechanical  contributions  from  vibrational  motion  to  energy  and  constant  pressure  heat  capacity. 
Theee  quantum  contributions  minus  the  classical  values  give  their  quantum  corrections.  They 
discuss  the  shift  in  the  vibrational  frequency  of  water  as  it  enters  the  liquid  phase  which 
changes  tha  zero  point  energy.  This  is  necessary  because  they  use  rigid  molecules.  The  type  of 
nonrigid  potential  which  we  use  includes  both  imra-  and  intermoiecular  degrees  of  freedom  and 
thus  in  principle  (but  not  yet  in  practice  due  to  potential  energy  function  inaccuracies  as  is  dis¬ 
cussed  below)  can  take  into  account  the  frequency  changes  from  gas  to  liquid. 

The  quantum  correction  technique  used  in  the  present  paper  involves  calculating  the  velo¬ 
city  spectrum  5(a)  from  molecular  dynamics  and  then  integrating  5(a)  over  all  frequencies 
with  a  weighting  function  which  is  the  difference  between  the  quantum  and  classical  harmonic 
weighting  functions  for  the  thermodynamic  variable  of  interest. 

A.  Velocity  spectrum 

The  velocity  spectrum  5(a)  of  a  classical  system  of  IV  atoms  in  equilibrium  is  defined  as 


in  which  v  is  frequency,  j8  -  (kg  T)~l  in  which  kg  is.  Boitzman's  constant  and  T  is  the  tempera¬ 
ture,  rnj  is  the  mass  of  the  atom  corresponding  to  the  jth  Cartesian  velocity  component  as  a 
function  of  time  vy(r),  and  <  >  indicates  an  average  over  the  ensemble.  The  spectral  density 
operator  D  (for  which  windowing  and  window  correction  techniques  are  described  eise- 
where-3*-  39 )  ^  evaluated  in  terms  of  probability  per  unit  angular  frequency, 

12 


•  f 

2>(v;(r)]  s  (2ir)~1lim~  jdt  exp(-i2irvt)vj(.t) 


(3.2) 


The  velocity  spectrum  may  also  be  computed  from  the  Fourier  transform  of  the  velocity  auto¬ 
correlation  function.  Note  that  the  velocity  spectrum  can  be  computed  separately  for  different 
subsets  of  atoms  (for  example,  different  elements,  different  chemical  environments  of  the  same 
element,  or  different  molecules)  and  the  velocity  spectrum  SM  can  then  be  computed  as  a 
sum  of  the  effects  from  these  different  subsets  of  atoms.  Thus,  as  we  will  see,  the  quantum 
corrections  also  can  be  partitioned  among  the  different  subsets  of  atoms.  Even  though  once  the 
dynamics,  i.e.  the  set  of  velocities  |vy(r)},  is  determined,  the  quantum  corrections  may  be  com¬ 
puted  separately  for  different  subsets  of  atoms,  it  should  be  remembered  that  normally  all 
atoms  together  contribute  to  determining  the  dynamics. 

m 

h  will  be  useful  below  to  know  the  value  of  the  integral  ^dvSfv).  The  Fourier 

transform  of  a  real  function,  e.g.  Vi(r),  has  an  even  real  pan  and  an  odd  imaginary  pan.40  The 
square  of  the  absolute  value  of  such  a  Fourier  transform,  e.g.  2>(vy(r)],  is  a  real  even  function. 
A  linear  combination  of  real  even  functions,  e.g.  S(v),  is  also  a  real  even  function.  Therefore 
S(-*)  -  S(v)  which  allows  us  to  write 

|dvSW  -  Jd*S(v)/2.  (3.3) 


Substituting  Eq.  (3.2)  into  Eq.  (3.1)  and  inserting  the  result  into  the  right  side  of  Eq.  (3.3) 
gives 


CdvSiv)  -  0  C  dv  mj<  lim  TT-  fdr  cxp(~l2wt)vj(t) 

*0  Zrn  Jm  |  *— m*T  St 


(3.4) 


Let 


»/  \  -T  <  /  <  r 

VJ(,) "  |0  otherwise, 

and  let  the  Fourier  transform  of  v/(r)  be  F}(r),  i.e., 

FJM  -  J*  txp(-/2rvi)  rf(t)  -  fdt  exp(-/2irw)  »y(r). 


(3.5) 

(3.6) 


Substituting  Eq.  (3.6)  into  Eq.  (3.4)  gives 

JdvS(v)  -  •  °-7) 

Exchanging  integration  and  the  r— ••  limit  gives 

<fw|/7(w)|*>  . 

By  Parse vai’s  theorem,40 


(3.8) 


•  6  • 


since  Vj(t)  is  a  real  function.  Substituting  Eq.  (3.9)  into  Eq.  (3.8)  gives 


J^dvSM  -  dt  |vy(r)j2> 

(3.10) 

■2L|vy(r)]2>  . 

Jm\ 

(3.11) 

Since  by  classical  equipartition  of  energy1 

(3.12) 

in  which  <£*>  is  the  average  kinetic  energy. 

f  dvS(v)  -  2ff(3JV/2ff)  -  3 N, 

(3.13) 

i.e.  the  integral  of  the  velocity  spectrum  from  zero  to  infinite  frequency  is  just  three  times  the 
number  of  atoms,  which  will  be  true  for  any  potential,  harmonic  or  not.  This  relation  can  be 
used  as  a  check  on  the  accuracy  of  computation  of  5(v),  and  to  interprete  S(.v)  in  terms  of  an 
equivalent  density  of  normal  modes  even  for  enharmonic  systems. 

The  diffusion  coefficient  has  a  particularly  simple  expression  in  terms  of  the  velocity  spec¬ 
trum.  The  diffusion  coefficient  3  of  a  particle  with  position  history  r(r)  is  defined  as41 

D  -  |  lim  ~  <  frW  “  r(0)l2>  (3.14) 

where  <  >  indicates  an  ensemble  average.  Letting  the  three  Cartesian  components  of  r (r)  be 
x(t),  y(t)  and  z(r),  we  have 

D  -  |  lim  < lx(r)  -  x(0)JJ  +  [y(r)  -  y(0)]J  +  U(r)  -  z(0)]2>.  (3.15) 

For  isotropic  systems,  the  equation  may  be  simplified  to 

B  -  lim  <  Ix(t)  -  jc(0)]j>,  (3.16) 

t — «•  2r 

or 

25  -  lim  <U(t)  -  x(— r)lJ>,  (3.17) 

where  x(r)  now  represents  any  one  of  the  three  Cartesian  components  of  r(r).  If  we  let 
D'lvjU)]  denote  the  value  of  the  spectral  density  at  zero  frequency,  then  Eq.  (3.2)  becomes 

D.[yj(t))  -  (2w)->  lim  ^-U(r)  -  x(-r)]2.  (3.18) 

Combining  Eqs.  (3.17)  and  0.18)  we  get 

25-w<$,[v>(/))>.  (3.19) 

V  SM  is  restricted  to  equivalent  particles,  then  Eq.  (3.1)  becomes 

5(e)-  </)tvy(r)]>  -  12wAfmff<D(vy(r))>  (3.20) 

.  y-i 

where  M  particles  are  being  considered  each  of  mass  m.  Then 

<$>,(*))>  -  S(0)f\2wMmfi  (3.21) 

and  thus  the  diffusion  constant  3  is  related  to  the  zero  frequency  value  of  the  velocity  spec¬ 
trum  5(0)  by 


-7- 


N 


5  -  S(0)/l2Mmfi,  (3.22) 

in  which  M  is  the  number  of  equivalent  particles,  m  is  their  mass,  and  0  -  (kB  T)~l  in  which 
kg  is  Boltzman’s  constant  and  T  is  the  temperature.  The  most  usual  application  of  Eq.  (3.22) 
is  to  consider  the  particles  to  be  molecules  and  to  compute  the  diffusion  constant  from  the  zero 
frequency  value  of  the  velocity  spectrum  of  the  center  of  mass  of  the  molecules. 


B.  Harmonic  approximation 

We  quantum  correct  the  classical  thermodynamic  variables  using  a  harmonic  oscillator 
approximation.  This  correction  is  based  on  a  division  of  the  dynamics  in  frequency  space.  The 
low  frequency  region  is  viewed  as  nearly  classical  but  containing  the  major  anharmonic  effects, 
and  the  high  frequency  region  is  viewed  as  nearly  harmonic  and  thus  can  be  quantum  corrected 
exactly  within  the  limits  of  the  harmonic  approximation. 

Consider  a  system  of  N  atoms  as  linked  by  harmonic  potentials. 

K<r*)  -  K0  +  i  £  Ar,Ar*  (3.23) 

in  which  A/j  and  A/*  are  displacements  from  a  potential  minimum  and  V0  is  the  potential 
energy  at  that  minimum.  Such  a  harmonic  situation  can  be  approached  classically  in  the  limit 
of  small  atomic  motions  about  a  potential  minimum,  i.e.  at  low  temperatures,  but  one  should 
remember  that  quantum  wave  functions  ample  the  potential  in  a  region  about  the  minimum 
even  at  absolute  zero,  and  thus  anharmonidty,  both  explicit  and  due  to  coupling  by  finite  dis¬ 
placements,  will  always  play  a  role  in  real  systems.  Nonetheless,  we  believe  that  at  higher  fre¬ 
quencies  an  analysis  which  uses  the  finite  temperature  classical  velocity  spectrum  interpreted  as 
if  it  were  fully  harmonic  will  usually  sufficiently  well  represent  the  thermodynamic  quantum 
corrections. 

In  the  harmonic  limit,  a  normal  mode  analysis  allows  us  to  view  the  system  as  a  set  of  3N 
harmonic  oscillators  with  qf  as  a  single  oscillator  partition  function.  The  total  canonical  parti¬ 
tion  function  Q  for  the  system  can  then  be  expressed  in  terms  of  the  partition  functions  qj  for 
the  individual  mottos  as 

(3.24) 

y-i 

or 

ln<?  -  £  tag).  (3.25) 

y-i 

If  the  normal  frequencies  are  continuously  distributed  we  may  take  the  integral 


tfoSGOln  q(p) 


(3.26) 


where  SM  is  the  density  of  normal  modes  with  frequency  v. 

To  show  that  the  velocity  spectrum  of  a  system  of  particles  linked  through  harmonic 
potentials  represents  the  density  of  normal  modes,  the  3  N  time  varying  Cartesian  position  com¬ 
ponents,  X], . . .  ,xj |V,  are  first  represented  in  terms  of  normal  coordinates.  We  have42 


qj  -  AjSin(ujt  +  9j)  (3.28) 

where  . 4j,v  are  the  normal  coordinates,  . «j.v  are  the  characteristic  normal 

mode  angular  frequencies  in  which  2 irvj  -  At  is  the  jth  normal  mode  amplitude.  9,  is  its 
phase,  and  a,k  are  constants  scaled  such  that 


reduced  energy  of  the  harmonic  oscillator,  0  -  (kg  T)~l,  h  is  Planck's  constant  and  v  is  the  fre¬ 
quency  of  the  oscillator.  Substituting  Eq.  (3.41)  into  Eq.  (3.26),  and  inserting  this  result  into 
Eqs.  (2.1)  through  (2.4)  gives 


-  Vo  +  kgT^dvSb)  Wgivk  wj(v)  -  1 

(3.42) 

m 

Cf  -  kg^ dvS(v)  (*);  (v)  -  1 

(3.43) 

m 

AC~V0  +  kgT^dvSiv)  Wfrh  WCA(V)  -  Inu 

(3.44) 

am 

Sc  ~  kg f  dvSiv )  Wf( v)\  W$(v)  -  [1  -  Inu]. 

il 

(3.45) 

These  classical  weighting  functions  Wc(v)  are  shown  in  Ftg.  1.  To  allow  the  zero  of  energy  to 
be  set  arbitrarily,  we  include  Vg,  the  energy  of  the  system  treated  classically  if  all  oscillations 
are  stilled.  The  expressions  for  energy  and  heat  capacity  reduce  to  the  familiar  classical  results 

P-Vg  +  lNkgT  (3.46) 

Cf  -  3 Nkg.  ,  (3.47) 


D.  Quantum  weighting  functions 

The  quantum  mechanical  partition  function  for  a  single  harmonic  oscillator  is  1 


1  -e- 


(3.48) 


where  the  superscript  Q  indicates  that  the  variable  is  derived  quantum  mechanically  and  again 
u  m  fihv  is  the  reduced  energy.  Substituting  Eq.  (3.48)  into  Eq.  (3.26)  and  inserting  this 
result  into  Eqs.  (2.1)  through  (2.4)  gives 


-  K0  +  kt  T^dvSir)  W°{v)\  W°(v)  -  |y  + 

C«  -  *a|*rf(v)  IK?>); 

A  0  -  V0  +  kgT^dvS(v)  WfWl  tvf(v)  -  | 

S  0  -  kt | dvS(v)  »f(«0;  tvfM  -  I ~  1 


-~r  -  In  (l  -  e-) 
e *  -  1 


(3.49) 

(3.50) 

(3.51) 

(3.52) 


Fig.  1  shows  these  quantum  weighting  functions  WQ(v). 

For  a  system  which  closely  approximates  a  set  of  harmonic  oscillators,  such  as  a  perfect 
crystal  at  low  temperature,1*9  the  above  equations  alone  can  be  used  to  compute  the  thermo¬ 
dynamic  variables. 

E.  Quantum  correction  weighting  functions 

The  quantum  corrections  (indicated  by  the  superscript  A  )  are  obtained  by  subtracting  the 
classical  representations  from  the  quantum  mechanical  representations  for  the  given  thermo¬ 
dynamic  variable. 

tr4(*)  .  iv°M  -  Wc(v) 


(3.53) 
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FIG.  1.  Universal  harmonic 
weighting  functions  W{v)  for 
the  thermodynamic  functions 
of  any  harmonic  system.  Dot* 
ted  lines  are  classical  Wc( v) 
from  Eqs.  (3.42)  to  (3.45), 
dashed  lines  are  quantum 
WQ{p),  from  Eqs.  (3.49)  to 
(3.52),  and  solid  lines  are 
quantum  correction 

Wx(v)  -  WQ(p)  -  Wc(r), 
from  Eqs.  (3.54)  to  (3.57).  It 
is  thus  the  solid  line  curves 
which  are  used  to  weight  the 
velocity  spectra  5(v)  to  com¬ 
pute  the  quantum  corrections 
to  the  thermodynamic  func¬ 
tions.  The  top  panel  is  for 
energy  £,  the  next  to  the  top 
panel  is  for  constant  volume 
heat  capacity  C„  the  next  to 
the  bottom  panel  is  for 
Helmholtz  free  energy  A,  and 
the  bottom  panel  is  for  entropy 
S.  The  lower  horizontal  scale 
is  plotted  in  reduced  energy 
u  -  j8 hv  in  which  £  -  kg  T,  kg 
being  Boltzman's  constant  and 
T  the  temperature,  while  h  is 
Planck’s  constant  and  v  the  fre¬ 
quency  of  the  oscillator.  The 
upper  horizontal  scale  is  the 
wavenumber  equivalent  to  u  at 
300  KL  Note  that  all  the  quan¬ 
tum  correction  weighting  func¬ 
tions  go  to  zero  at  low  fre¬ 
quency  where  enharmonic 
effects  become  important. 
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£A .  £9  -  £c  -  k§  rjdt'Sh)  Wg(vY, 

fV?(v)  - 

2  e"— 1 

k*|dvS(v)  W^(v); 

Kiv)  “ 

[ 

l  (1  -  e“)J 

m 

A*-AQ-AC-  kBT^dvS(v) 

*rfw  - 

1  -  e- 
ln  e-un  ~ 

(3.54) 

(3.55) 

(3.56) 


S4  -  SO  -  S*  -  kB] dvSiv)  J V*(v)\  W/(*)  - 


eu-  1 


-  Ml  -  e~“) 


+  Inu  -  l).  (3.57) 

The  quantum  correction  weighting  functions  W*(v)  are  also  shown  in  Fig.  1.  Note  that  follow¬ 
ing  Eqs.  (3.1),  (3.2),  and  (3.54)-(3.57)  we  can  partition  if  we  wish  the  quantum  corrections 
among  different  subsets  of  atoms,  for  example  different  elements,  different  chemical  environ¬ 
ments  of  the  same  dement,  different  molecules,  or  molecules  in  different  environments. 


IV.  MOLECULAR  DYNAMICS  RESULTS  FOR  LIQUID  WATER 

Water  is  the  most  important  of  all  solvents,  and  the  molecular  level  understanding  of  its 
bulk  properties  is  of  considerable  intrinsic  interest.  We  have  thus  chosen  it  as  a  test  case  for 
our  techniques.  A  quantum  calculation  for  a  system  of  molecules  large  enough  to  adequately 
represent  liquid  water  is  at  present  impractical,  and  thus  thermodynamic  quantities  are  com¬ 
puted  by  classical  mechanics,  usually  by  Monte  Carlo  or  molecular  dynamics  techniques.  Such 
classical  molecular  mechanics  calculations  on  liquid  water  have  been  discussed  in  reviews  by 
Stillinger,35  Barnes,44  Wood,45  and  Beveridge  et  al.46  God  and  Hockney47  have  written  a 
comprehensive  bibliography  for  earlier  molecular  dynamics  in  general.  It  win  be  shown  for 
liquid  water  that  quantum  corrections  are  needed  for  both  inter  and  intramolecular  motions  to 
match  experimental  quantum  reality. 


A.  Liquid  water  potentials  and  previous  computer  sfanulatloas 

A  major  obstacle  for  any  molecular  mechanics  computer  simulation  is  the  development  of 
an  accurate  potential  surface.  Both  experimental  data  and  quantum  calculations  are  valuable  to 
tins  end.  Bernal  and  Fowler4*  (BF)  in  1933  have  given  a  rigid  three  point  charge  plus 
Leonard- Jones  potential  for  water.  An  empirical  potential  for  water  was  introduced  in  19S1  by 
RowHnson49  (ROW),  and  tested  by  Barker  and  Watts.5®-52  This  is  a  rigid  four  point  charge 
model  with  a  Lennard-Jones  oxygen-oxygen  potential.  The  analytical  form  of  the  Rowlinson 
potential  has  been  utilized  in  several  improved  potentials,  namely  BNS  and  ST2.  Ben-Naim  and 
Stillinger  introduced  the  BNS  potential55  in  1972,  and  Rahman  and  Stillinger54*5*  and  others52 
utilized  it  in  several  test  studies.  After  finding  the  potential  too  tetrahedronally  directional55 
and  noting  an  improvement  after  an  energy  rescaling54  Stillinger  and  Rahman  introduced  die 
ST2  potential57  in  1974.  ST2  has  been  used  and  tested  extensively  by  many  workers.5*47 

Shipman  and  Scheraga**  have  developed  a  seven  point  charge  effective  pair  potential  (SS) 
for  water  using  a  variety  of  experimental  data.  Both  water  dimer**  and  ice-like  water  duster70 
studies  have  been  carried  out.  In  an  attempt  to  indude  intramolecular  vibrations,  to  allow  for 
possible  molecular  dissociation,5*  and  to  account  for  some  of  the  nonadditive  interaction71 
between  waters,  Lemberg  and  Stillinger  introduced  a  central  force  potential19  (LS)  in  1975.  In 
this  scheme  both  bonded  and  non  bonded  oxygen-hydrogen  interactions  use  the  same  potential 
as  do  all  hydrogen-hydrogen  interactions.  It  has  been  both  ftmher  applied72  and  improved75 
(LS2). 

An  ab  Inttto  water  potential  prepared  by  fitting  analytical  functions  to  quantum  mechani¬ 
cally  calculated  energy  versus  nudear  position  data  was  developed  by  Popkie.  Kistenmacher  and 
dementi74  in  1973  at  the  Hartree-Fock  level.  Several  studies75*80  have  used  this  potential 


(HF),  and  ithas  been  found  that  the  neglect  of  electron  correlation  effects  leads  to  significant 
inaccuracies.'5*77'7*  In  response,  Matsuoka,  Clementi  and  Yoshimine*1  carried  (Hit  ab  initio  cal¬ 
cinations  with  configuration  interaction  to  account  for  correlation  effects.  The  resulting  poten¬ 
tial  (Cl)  has  been  tested  by  several  investigators.65*67**2*87  A  specific  comparison  of  the  HF  and 
Cl  potentials  has  been  carried  out  by  Swaminathan  and  Beveridge.**  The  Cl  potential  has  been 
criticized  because  of  lack  of  agreement  with  experimental  liquid  density,57  its  poor  reproduction 
of  the  second  viriai  coefficient  of  steam*5  and  the  high  rms  error  in  fitting  the  original  calcula¬ 
tions.*9 

Watts  has  provided  a  flexible  water  dimer  potential90  (WATTS)  in  which  a  largely  empiri¬ 
cal  imermolecular  potential  is  complemented  by  an  intramolecular  potential  derived  from  vibra¬ 
tional  spectroscopy.  The  WATTS  potential  has  been  studied  and  criticized  by  McDonald  and 
Klein.  *9*91 

Several  more  recent  water  potentials  deserve  mention.  Stiitinger  and  David92  have 
developed  a  polarization  potential  (SD)  in  order  to  describe  deformable  water  molecules  and 
Stiitinger95  has  studied  its  dynamic  properties.  The  importance  of  non  pairwise  additive  effects 
for  water  is  also  stressed  by  Barnes  et  al,94  and  they  introduce  a  polarizable  electropole  model 
(PE)  and  test  it  in  a  water-amino  add  system.95  Goodfellow96  continues  this  study  of  coopera¬ 
tive  effects  and  the  PE  potential  by  showing  how  the  PE  model  can  be  effldently  applied.  The 
present  parameterization  of  the  PE  model  has  been  criticized  recently  because  the  oxygen- 
oxygen  radial  distribution  agrees  poorly  with  experiment46  Nemenoff,  Snir,  and  Sdieraga97 
have  developed  an  empirical  technique  (EPEN/2)  for  potential  function  development  which  has 
been  revised  by  Marchese,  Mehrotra  and  Beveridge.9*  Berendsen  et  al99  have  produced  a  single 
point  charge  (SPC)  potential  with  Lennard-Jones  interaction  between  oxygens  in  order  to  han¬ 
dle  conveniently  protein-water  systems.  Jorgensen  has  developed  a  set  of  transferable  inter- 
molecular  potential  functions  (TTPS)  for  application  to  organic  liquids  and  water.100  Reimers, 
Watts,  and  Klein101  propose  a  revised  version  (RWK2)  of  the  WATTS  potential. 

Many  molecular  dynamic54* 57*59^1*64*72*86**7*91* 102  calculations  have  been  carried  out  as 
well  as  Monte  Carlo7,  **57, 461 5^52*62*65,67*74*75*7***2**4**5*100*101* 105,104  ^jniiitiw**  on  liquid 
water  using  most  of  the  potentials  described  above.  In  addition,  Weres  and  Rice105  discuss  the 
calculation  of  liquid  water  thermodynamic  properties  and  their  quantum  corrections  using  the 
BNS  potential  (with  some  modifications)  and  a  cell  model  viewpoint. 

Several  papers  have  tested  and  compared  the  variety  of  water-water  potentials,  often  with 
disappointing  results.  Morse  and  Rice106  calculate  some  of  the  properties  of  ice  with  many  of 
the  above  potentials.  The  results  raise  serious  questions  about  the  ability  of  ST2,  WATTS,  and 
LS2  to  accurately  reproduce  the  properties  of  ice  while  Cl,  with  the  exception  of  reproducing 
too  low  densities,  fates  wail.  Reimers  and  Watts101  make  a  related  comparison  extending  to  all 
three  phases.  WATTS,  ROW,  and  BNS  reproduce  the  second  viriai  coefficient  of  steam  well 
but  fail  in  the  other  two  phases.  Cl  and  ST2  do  well  for  the  liquid  phase  but  fail  with  ice  and 
steam.  They  conclude  that  all  models  tested  are  generally  inadequate  u>  handle  all  three 
phases,  but  that  their  revised  RWK2  potential  fares  best 

B.  Molecular  dynamics 

Our  molecular  dynamics  calculations  are  carried  out  on  a  system  of  250  water  molecules  at 
a  density  of  1.0  g  cm*5  and  a  temperature  of  300  K  with  cubic  periodic  boundary  conditions 
using  a  special  molecular  mechanics  package  running  on  an  array  processor.107- 108  Experimen¬ 
tally,  this  density  corresponds  to  a  pressure  of  85  atm  with  a  negligible  resulting  difference109  of 
0.012  kJ  mole*1  in  total  energy  compered  to  a  pressure  of  1  atm  which  corresponds  to  a  density 
of  0.997  g  cm*5.  Previous  molecular  dynamics  calculations  of  thermodynamic  quantities  for 
water  have  been  carried  out  using  an  array  processor  by  Rapapon  and  Sdieraga17* 1,0  who  stu¬ 
died  a  sample  of  343  rigid  waters  using  the  G  potential  with  long  runs  and  by  Swope,  Ander¬ 
sen,  Berens,  and  Wilson102  who  studied  the  properties  of  water  clusters.  The  software  used 
previously107* 108  has  been  augmented  by  an  imermolecular  force  and  energy  calculator  for 
water  as  implemented  by  Swope  and  Andersen.111  This  calculator  utilizes  a  piecewise  fifth  order 
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FIG.  2.  Velocity  spectra  times  the  speed  of  light  c  normalized  for  one  molecule  of  HjO  at  300 
K  and  1.0 1  cm-1,  using  the  Wans  potential  with  2S0  waters  and  cubic  periodic  boundary  con¬ 
ditions.  The  lower  panel  contains  c S*(v)  for  the  hydrogen  atoms  per  molecule  of  water, 
eSo(v)  for  the  oxygen  atoms  per  molecule  of  water,  the  sum,  cSifM  +  cSoM  •  and 

the  center  of  mass  velocity  spectrum  cScm(v).  The  upper  panel  is  a  blowup  of  the  low  fre¬ 
quency  region  of  the  lower  panel.  The  lower  horizontal  scales  are  in  terms  of  the  reduced 
energy  u  -  &hv  in  which  0  -  (kg  D'1  and  v  is  the  frequency.  The  speed  of  light  c  is  included 
so  that  the  integral  of  eS(*»)  in  cm  vs  the  upper  scale  of  wavenumbers  in  cm*1  will  be  dimen¬ 
sionless,  giving  the  total  number  of  equivalent  harmonic  oscillators.  For  a  purely  harmonic  sys¬ 
tem  the  velocity  spectrum  5(»)  would  give  the  number  of  harmonic  modes  per  unit  frequency. 
Note  that  the  H  atoms  dominate  Sn} o(»)  above  300  cm'1  and  the  O  atoms  below  it. 


polynomial  fitted  to  the  analytical  potential  energy  functions  as  a  function  of  the  square  of  the 
distance  between  the  two  atoms  being  considered.  It  thus  both  allows  a  general  algorithm  to 
evaluate  the  polynomial  previously  fitted  to  arbitrary  analytic  functions  and  eliminates  the 
necessity  of  a  square  root  operation. 

The  method  for  applying  a  switching  function  as  developed  by  Andersen  and  Swope 
smooths  each  water-water  energy  contribution  to  zero  as  the  corresponding  oxygen-oxygen  dis¬ 
tance  passes  through  the  switching  region,  which  for  our  system  extended  from  0.85  to  0.90 
nm.  This  technique  eliminates  the  problem  of  artificially  created  monopoies  (and  possibly  large 
dipoles)  normally  encountered  by  an  atom-atom  force  feathering  or  truncation  method  as  only 
pan  of  the  water  molecule  passes  through  the  feathering  region  (and  is  possibly  imaged).  This 
artifact  is  especially  pronounced  with  water  unless  the  Andersen-Swope  technique  is  used,  as 
the  partial  charges  on  each  atom  are  relatively  large. 

The  semiempirical  flexible  water  molecule  potential  developed  by  Witts90  is  used.  The 
intramolecular  potential  is  a  standard  Taylor’s  series  in  internal  coordinates  about  the  potential 
minimum  as  derived  from  vibrational  spectroscopy.  The  intertnoiecular  potential  is  pairwise  by 
atoms  ami  fitted  to  the  second  virial  coefficient  of  steam. 

Equilibration  of  die  initial  water  system  is  accomplished  by  following  periods  of  dynamics 
(0.1  •  10  picoseconds)  with  randomizations  of  velocity  according  to  a  MaxweU-Boltzmann  dis¬ 
tribution  at  the  desired  temperature  until  the  temperature  of  the  system  stabilizes.  The  total 
simulation  time  involved  in  equilibration  is  approximately  60  picoseconds.  The  time  step  of 
integration  during  equilibration  is  0.5  femptoseconds  while  for  the  data  collection  a  time  step  of 
0.25  femptoseconds  is  used. 

The  velocity  data  is  accumulated  by  selecting  out  the  velocities  every  12  time  steps  over  a 
period  of  50  000. time  steps  (113  picoseconds  total  simulation  time).  A  more  elegant  approach 
would  be  to  use  a  digital  tow-pass  filter  before  sampling.113  The  energy  and  heat  capacity  data 
are  the  result  of  a  much  longer  series  of  seven  runs  for  a  total  of  380  000  time  steps  over  95 

C.  Velocity  spectra 

The  velocity  spectra  S(v)  shown  in  Fig.  2  are  calculated  by  fast  Fourier  transforms  of  the 
velocity  time  histories  of  various  components  of  the  system.  We  define  the  following  normal¬ 


ized  velocity  spectra 

S„M «/'<Z){v/'(/)]>  (4.1) 

M  y- 1 

SoM  -  £  mf><D[vf>U))  >  (4.2) 

M  y-i 

SffjO  W  •  5p(r)  +  SffM  (4.3) 

Sarto  -  n-f  £  \>  (4.4) 

A*  Jm  | 


where  m°,  m M,  and  mHi°  represent  the  masses  of  an  rnogen  atom,  hydrogen  atom,  and  water 
molecule  respectively;  D  is  the  spectral  density  operator  defined  in  Eq.  (3.2);  v/fy),  and 
vP*(r),  represent  the  velocity  time  histories  of  the  jth  oxygen  atom,  hydrogen  atom  and  molec¬ 
ular  center  of  mass  respectively;  Af  is  the  number  of  water  molecules,  where  Af  s  W3;  and  a 
factor  of  1  /Af  has  been  introduced  to  normalize  the  velocity  spears  to  that  for  one  molecule  of 
water.  The  contribution  to  by  both  the  oxygens  and  the  hydrogens  is  determined  by 

computing  each  spectrum,  5o(*)  and  5/y(»),  separately.  The  high  frequency  vibrational  peaks 
composed  mainly  of  the  oxygen-hydrogen  vibrations  are  easily  seen  in  Fig.  2.  The  center  of 
mass  velocity  spearum  of  the  system  is  also  computed  and  its  spectrum  reflects  the  highly 
damped  vibrational  modes  of  whole  water  molecules. 


The  area  under  Sh:oM  in  Fig.  2  equals  9.0,  the  equivalent  number  of  harmonic  oscilla¬ 
tors  (including  hindered  translation  and  rotation)  per  molecule  of  water,  as  expected  from  Eq. 
(3.13).  (The  speed  of  light  is  introduced  to  make  the  integral  versus  cm-1  unitless.)  The  dou¬ 
ble  peak  in  the  range  2600-5000  cm**1  which  corresponds  to  the  symmetric  and  asymmetric 
stretching  modes  of  the  water  molecule  has  an  area  of  1.S9.  The  peak  in  the  range 
1200-2600  cm'1  which  corresponds  to  the  bending  of  the  HOH  bond  angle  has  an  area  of  1.00. 
This  substantiates  the  view  of  SM  as  a  density  of  normal  modes  and  further  suggests  that  the 
dose  association  of  the  water  molecules  in  the  liquid  phase  h  \s  shifted  some  of  the  high  fre¬ 
quency  stretching  motion  down  into  the  low  frequency  region.  * 

In  prindpie  a  potential  with  both  intra  and  intermoiecular  degrees  of  freedom  like  the 
WATTS  potential  we  have  used  coukl  take  into  account  the  frequency  changes  from  gas  to 
liquid.  The  actual  frequencies  for  the  WATTS  potential  for  the  gas  phase  should  be  dose  to 
the  harmonic  values11*  of  »\  •  3132  cm*1  (symmetric  stretch),  v2  -  1649  cm'1  (bend),  and 
-  3943  cm'1  (asymmetric  stretch),  compared  to  the  computed  liquid  phase  peaks  centered  at 
3680,  1740,  and  3760  cm'1  as  shown  in  Fig.  1  In  real  water,  the  infrared  and  Raman  spectra 
show  the  gas  phase  anharmonic  frequencies114  to  be  3652,  1595,  and  3756  cm'1  and  the  liquid 
phase10- 112  shows  a  bending  peak  at  approximately  1650  cm'1  and  a  brood  stretching  peak  cen¬ 
tered  at  approximately  3400  cm'1  with  perhaps  a  subsidiary  peak  at  approximately  3200  cm'1 
Thus  the  WATTS  vibrational  shifts  from  gas  to  liquid  phase  qualitatively  resemble  the  real 
water  shifts  with  large  shifts  downward  in  frequency  for  the  stretching  motions  and  a  smaller 
shift  upward  for  the  bending  motion,  but  the  agreement  is  certainly  not  quantitative. 

From  5o*(0)  in  Fig.  2  and  Eq.  (3.22)  we  obtain  for  the  censer  of  mass  diffusion 
coefficient  3  of  water  a  value  of  4.01  x  I0r*mh~l  compared  to  the  experimental  value116- 117  of 
142  x  10-,mV1  for  liquid  water  at  300  KL  The  precision  of  our  reported  value  is  questionable 
because  we  selected  out  every  twelfth  velocity  rather  than  all  velocities  for  the  fast  Fourier 
transform  due  to  computer  memory  Ihnharions,  and  a  more  reliable  value  could  be  computed 
from  the  asymptotic  slope  of  the  mean  square  displacement  of  the  center  of  mass  for  a  long 
molecular  dynamics  run.  It  should  also  be  remembered  that  the  finite  size  of  the  periodic 
boundaries  may  aflbet  the  longest  wavelength  and  lowest  frequency  motions  and  in  particular 
that  hydrodynamic  or  concanad  motions  involving  many  moiacuias  may  not  be  accurately  han¬ 
dled. 

Berendsen  et  el*  have  reportad  a  spectral  density  of  the  center  of  mass  of  rigid  molecule 
liquid  water,  using  the  previously  dmcrihed  SPC  potential,  which  is  strikingly  similar  to  our 
SatM.  They  report  a  diffusion  coefficient  of  3.6  x  10"Ws'1. 

D.  Quantum  censcdsne 

The  difference  between  the  classical  and  quantum  mechanical  weighting  functions  W(v) 
arises  from  the  difference  between  the  classical  and  quantum  harmonic  oscillator  partition  func¬ 
tions  q(v).  hr  the  damtari  limit  of  ft  —  0,  or  equivalently  »-0,r-0,or  T  —  «•,  this  dis¬ 
tinction  disappears, 

Um  q°(v)  -  Um  qc(v).  (4.5) 

*  — 0  *  — 0 

This  implies 

Um  WQ{v)  -  Um  WQ(V)  -  Um  WCM  (4.6) 

A  — 0  »  — 0  »  —  0 

and  thus, 

1Fa(0)  -  WQ(0)  -  IVC(0)  -  0  (4.7) 

in  all  cases,  as  can  be  seen  in  Fig.  1.  The  divergence  of  WQ{ »)  from  Wc(v)  as  v  increases 
results  in  a  preferential  weighting  of  high  frequency  motions  in  the  calculation  of  quantum 
corrections. 


Table  I  gives  the  liquid  water  quantum  corrections  computed  from  the  velocity  spectrum 
5w,o(»»)  from  Eqs.  (4.1)  to  (4.3)  as  shown  in  Fig.  2  and  the  quantum  correction  weighting 
functions  W*(»)  as  shown  in  Fig.  1,  using  Eqs.  (3.54)  to  (3.57).  The  curves  of  the  products  of 
S(v)  and  W*  for  energy,  heat  capacity,  free  energy  and  entropy  are  shown  in  Figs.  3  to  6, 
illustrating  the  contribution  to  the  quantum  corrections  as  a  function  of  frequency  and  of  atom 
type.  A  separation  is  made  for  the  purposes  of  Table  I  in  frequency  space  at  1200  cm*1 
between  the  intermotecuiar  and  intramolecular  motions  for  liquid  water.  Norn  that  the  inter* 
molecular  motion,  the  hindered  translation  and  rotation,  contribute  substantially  to  the  total 
quantum  corrections. 


TABLE  L  Inter-  and  intramolecular  contribution  to  liquid 
water  quantum  correction  at  300  K,  per  mole. 


Inter 

(0-1200  cm*1) 

Intra 

(1200-5000  cm*1) 

Total 

£A(kJ) 

4.2 

45.0 

49.2 

C*(J/K) 

•11.0 

•23.8 

•34.8 

/4AtkJ) 

2.2 

33.4 

35.6 

SaCJ/K) 

6.5 

318 

45.3 

Classically,  a  harmonic  oscillator  contributes  kgT  to  the  energy  regardless  of  frequency  as 
a  result  of  equipertition  of  energy.  This  produces  a  straight  line  for  the  classical  weighting 
Auction  in  the  top  panel  of  Fig.  1.  Quantum  mechanics,  however,  requires  that  a  harmonic 
oscillator  contain  a  minimum  or  pm  point  energy  of  hv/ 2.  For  a  harmonic  oscillator  with 
hv  «  kgT,  this  requirement  is  unimportant  and  quantum  effects  are  small.  In  contrast,  a 
quantum  harmonic  oscillator  with  hv  »  kgT  has  an  average  energy  near  h»/2.  As  a  result 

plto_  Wftv)  -w/2-g^L.  (4.8) 

Thus  the  quantum  effects  are  large  for  a  high  frequency  harmonic  oscillator  as  it  contributes 
hv/ 2  to  the  energy  instead  of  kgT.  Table  I  shows  a  value  of  49.2  kJ  for  the  total  quantum 
correction  to  energy.  Others  have  accounted  for  this  quantum  effect  by  introducing  a  constant 
into  the  potential  energy  Auction.  Using  spectroscopic  data,  Eieenberg  and  Kauzmann10  have 
calculated  S5.4S  kJ  as  a  zero  point  energy. 

Heat  capacity  is  unique  in  that  h  results  in  a  negative  quantum  correction,  and  it  has  the 
most  significant  contribution  front  the  low  frequency  region  compared  to  the  other  corrections 
we  have  listed.  As  a  remit  of  aqui  partition  of  energy,  the  classical  harmonic  oscillator  contri¬ 
butes  kg  to  C„  regardtau  of  frequency.  This  produces  a  straight  line  in  the  next  to  the  top 
panel  of  Fig.  1.  In  contrast,  the  qumtum  mechanical  harmonic  oscillator  with  hv  »  kgT  is 
"stuck"  in  the  ground  state  and  changes  very  little  in  response  to  changes  in  temperature.  As  a 
result, 

Um  W*  (v)  -0.  (4.9) 

» — «•  *■* 

Thus,  for  each  harmonic  oscillator  with  hv  »  kg  T.  kg  mm  be  subtracted  from  the  classically 
calculated  C«.  The  importance  of  the  low  frequency  contribution  to  the  quantum  correction  for 
constant  volume  heat  capacity  results  from  the  rapid  divergence  of  W^(v)  and  »£>>  “  ' 
increases  from  zero. 

The  equation  A  •  £  -  75  holds  in  an  analogous  manner  for  the  quantum  corrections  as 
a  result  of  the  linear  form  of  the  quantum  correction  equations.  The  energy  terra  dominates 
for  harmonic  oscillators  and  thus  the  quantum  correction  for  free  energy  is  always  positive. 

The  reader  may  be  surprised  that  the  quantum  correction  for  entropy  is  positive.  A  quan¬ 
tum  mechanical  harmonic  oscillator  with  hv  »  kgT  is  stuck  in  the  ground  state  and 
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.  4.  Constant  volume  heat  capacity  quantum  correction  curves  for  liquid  water  for  H  atoms, 
0  atoms,  and  their  sum  giving  the  total  HjO  quantum  correction.  Plotted  is  the  product  of 
speed  of  light  c,  the  velocity  spectrum  5(»),  and  the  heat  capacity  quantum  correction 
hting  function  (»»)  vs  the  reduced  oscillator  energy  umfihv  on  the  bottom  axis  and 
wavenumber  equivalent  at  300  K  on  the  top  axis.  The  integral  of  the  product  SM  {») 
vs  w  gives  the  quantum  correction  to  constant  volume  heat  capacity,  as  shown  in  Eq.  (3.55). 


FIG.  6.  Entropy  quantum  correction  curves  Tor  liquid  water  for  H  atoms,  for  O  atoms,  and  their 
sum  giving  the  total  HjO  quantum  correction.  Plotted  is  the  product  of  the  speed  of  light  c, 
the  velocity  spectrum  5(#>),  and  the  entropy  quantum  correction  weighting  function  W$(»)  vs 
the  reduced  oscillator  energy  uMfihv  on  the  bottom  axis  and  the  wavenumber  equivalent  at 
300  K  on  the  top  axis.  The  integral  of  the  product  SM  #'/(»»)  vs  v  gives  the  quantum  correc¬ 
tion  to  entropy,  as  shown  in  Eq.  (3.57). 


(4.10) 


contributes  almost  nothing  to  the  entropy.  Thus,  as  seen  in  Fig.  1, 

lim  -  0. 

*  —  «• 

to  contrast,  the  classical  harmonic  oscillator  weighting  function  has  the  following  properties. 

lim  wf(p)  -  oo.  (4.H) 

lim  WfM  -  (4.12) 

The  first  equation  indicates  that  an  unconstrained  particle  has  an  unlimited  number  of  available 
states.  The  second  equation  results  from  the  difficulty  of  applying  the  third  law  of  thermo* 
dynamics  to  the  classical  representation  of  entropy  for  a  harmonic  oscillator.  Because  of  the 
negative  sign  of  the  classical  weighting  function,  the  quantum  correction  for  entropy  is  positive. 

Figures  3  to  6  show  the  products  of  the  velocity  spectra  S( v)  with  the  quantum  correction 
weighting  functions  WA(v)  for  energy  E,  constant  volume  heat  capacity  C„,  Helmholtz  free 
energy  A,  and  entropy  S.  Since  SHto (*)  can  be  partitioned  into  separate  contributions  from 
the  hydrogen  and  oxygen  atoms,  we  also  partition  the  products  (v)  W1[p)  and  thus  com¬ 
pute  separately  the  hydrogen  and  oxygen  atom  contributions  to  the  quantum  corrections.  The 
hydrogen  atom  motions  dominate  except  at  the  very  lowest  frequencies  which  nave  little  weight 
anyway. 

E.  Energy 

Seven  water  samples  with  different  energies  are  created  and  equilibrated,  and  the  average 
temperature  for  each  sample,  calculated  over  at  least  ten  picoseconds  running  time,  is  plotted  in 
Fig.  7.  A  straight  line  is  fitted  to  the  points,  and  the  total  classical  energy  E0  corresponding  to 
300  K  is  calculated.  By  averaging  over  a  subset  (500  time  steps  selected  over  a  time  period  of 
1.25  picoseconds)  of  a  complete  run  at  300  K  we  also  compute  the  average  value  of  the 
intramolecular  potential  energy  K**.;  the  intermolecular  potential  energy  K**;  and  the  kinetic 
energy  &  as  shown  to  Table  II.  Because  Ec  is  calculated  from  the  fitting  shown  in  Fig.  7  while 
and  Ek  are  calculated  from  the  short  subset  discussed  above,  there  is  a 
0.1  kJ  mole*1  discrepancy  between  the  values  shown  in  Table  II  for  Ec  and  for  the  sum  of  its 
components  +  V^-b  Ek-  Die  quantum  correction  Ex  is  obtained  by  integrating  the 
function  SM  Wj(y)  as  shown  in  El  (3.54)  and  Fig.  3.  Addition  of  the  kinetic  energy  (calcu¬ 
lated  from  instantaneous  velocities10*11* )  to  the  total  potential  energy  results  in  conservation 
of  energy  to  one  part  in  thirty  thousand  with  the  0.25  fs  integration  step  size  used.  As  sug¬ 
gested  by  Andersen119  we  graphed  the  standard  deviation  of  the  total  energy  versus  the  square 
of  the  time  step  for  several  molecular  dynamics  runs.  The  resulting  nearly  linear  plot  verifies, 
the  accuracy  of  our  software  and  hardware  as  the  Veriet  integration  algorithm  gives  an  error  in 
total  energy  in  proportion  to  the  square  of  the  integration  time  step.119  To  calculate  K**,  we 
first  remove  a  constant  representing  the  zero  point  energy  contribution  from  the  original 
WATTS  intramolecular  potential.90 


TABLE  II.  Energy  (kJ  mole-1). 

Vm 

5.2 

Vmm 

-42.1 

Ek 

11.2 

& 

— 25.6a 

£A 

49.2 

E*~  -  &  +  £A 

23.6 

£■* 

21.5* 

'Calculated  from  Fig.  7 
‘See  Table  IV. 


McDonald  and  Klein19  calculate  with  molecular  dynamics  an  internal  potential  energy  of 
—33  kJ  mole'1  for  the  WATTS  potential  for  273  K  and  1  gm  cm-1  and  Reimers  and  Watts101 
report  a  Monte  Carlo  calculation  giving  an  internal  energy  of  -29.2  kJ  mole*1  for  the  Watts 
potential  at  298  K  and  0.997  gm  cm'1  Both  calculations  differ  from  the  present  one  in  that  their 
water  molecules  are  constrained  to  be  rigid.  To  determine  the  effect  of  this  we  increased  the 
force  constants  of  our  waters  first  by  a  factor  of  four  and  then  sixteen  while  decreasing  the  time 
step  first  by  two  and  then  by  four.  The  result  is  that  the  intermolecular  potential  energy 
decreases  (becomes  more  negative)  with  changes  on  the  order  of  1  kJ  to  2  kJ,  indicating  that 
the  introduction  of  flexible  waters  increases  (nukes  less  negative)  the  intermolecular  potential 
energy  over  a  rigid  water  calculation.  Reimers  and  Watts  run  at  298  K  and  the  1  atm  density  of 
0.997  gm  cm'1  compared  to  our  temperature  of  300  K  and  85  atm  density  of  1.0  gm  cm'1.  We 
performed  a  special  test  run  at  Q.977  gm  cm'1  and  298K  and  calculated  an  intermolecular 
potential  energy  only  marginally  different  from  the  1.0  gm  cm'1  value,  in  line  with  the 
0.012  kJ  mole'1  shift  expected199  for  the  total  potential  energy  from  experimental  thermo¬ 
dynamic  measurements.  We  perform  a  molecule-by-molecule  imaging  with  force  feathering 
technique  following  Andersen111  rather  than  potential  cut-offs  as  used  by  Reimers  and  Watts  or 
Ewaid  sums  as  used  by  McDonald  and  Klein.  To  explore  the  effects  of  potential  or  force 
smoothing  or  cutoff,  we  carried  out  several  additional  test  runs  whose  results  are  summarized 
in  Table  ID.  The  standard  deviations  are  given  within  parentheses  and  a  time  step  of  0.25  fs  is 
used  for  each  ran. 


TABLE  m.  Energies  in  kJ  mole'1  (with  standard  deviations 
given  in  parentheses),  as  well  as  bond  length  and  bond  angle  distortions 
_ for  several  cutoff  and  feathering  boundary  methods. 


T>pe 

y<mm 

y*m 

& 

~  (%) 

rm 

(%) 

»«* 

ANDERSEN 

5.2  (0.22) 

•411 

(0.26) 

-25.6 

(0.00099) 

0.52 

-1.0 

CUTOFF  I 

5.5  (0.32) 

•58.8 

(6.1) 

-41.8 

(6.4) 

0.56 

-1.0 

CUTOFF  n 

5.5  (0.32) 

•50.4 

(0.45) 

-33.6 

(0.029) 

0.56 

•1.0 

AA  SMOOTH 

27.5 

•111.7 

1.8 

-1.7 

IO 

4.7 

0.2 

-0.1 

Boundary  effects  are  a  significant  problem  for  systems  like  liquid  water  where  the  long 
range  Coulombic  forces  extend  well  beyond  the  dimensions  of  the  model.  One  way  to  deal 
with  these  nonzero  forces  near  the  boundary  is  to  choose  a  cutoff  distance  beyond  which  the 
potential  energy  is  set  to  zero.  For  an  atom-atom  central  force  system,  this  cutoff  of  the  poten¬ 
tial  can  also  be  done  atom  by  atom  (CUTOFF  I  in  Table  III),  and  the  resulting  forces  necessary 
for  molecular  dynamics  calculations  are  then  the  derivatives  of  the  potential  within  the  cutoff 
distance  and  zero  beyond,  with  a  delta  function  at  the  boundary  which  being  of  measure  zero  in 


length  is  never  seen  by  the  dynamics  calculation.  Such  energy-force  pairs  are  inconsistent  due 
to  the  effective  neglect  of  the  delta  function  force  term  at  the  cutoff  distance  which  prevents 
conservation  of  energy  in  actual  molecular  dynamics  runs  as  required  for  microcanonical  sys¬ 
tems.  The  atoms  fail  to  feel  the  force  delta  function  and  can  drift  back  and  forth  over  the 
cutoff  boundary  with  resulting  large  potential  energy  fluctuations.  For  ordinary  Monte  Carlo 
systems  where  forces  are  not  needed,  this  difficulty  is  avoided.  We  suspect,  however,  that  the 
large  fluctuations  in  the  radial  distribution  function  which  occur  at  the  cutoff  distance  introduce 
significant  perturbations  to  the  system.  Table  III  shows  the  results  of  a  sample  molecular 
dynamics  calculation  (CUTOFF  I)  using  a  cutoff  of  0.875  nm  at  the  midpoint  of  th  0.85  to 
0.90  nm  Andersen-Swope  feathering  which  we  used  for  our  actual  thermodynamic  calculations. 
Notice  that  the  standard  deviation  (in  parentheses)  for  the  intermolecular  potential  energy  is  a 
full  fifty-seven  percent  the  kinetic  energy. 

By  cutting  off  the  entire  potential  molecule-by-molecuie64  using  either  the  distance 
between  the  two  centers  of  mass  or  the  very  similar  oxygen-oxygen  distance  as  the  functional 
variable,  the  effectively  neglected  delta  function  force  terms  in  molecular  dynamics  calculations 
are  reduced  significantly  in  magnitude  as  they  now  represent  truncated  dipole-dipole  rather  than 
monopole-monopole  interactions.  Molecule-by-molecule  cutoff  is  also  preferable  to  atom-by¬ 
atom  cutoff  for  Monte  Carlo  calculations  as  the  fluctuations  in  the  radial  distribution  function  at 
the  cutoff  distance  should  be  greatly  reduced. 

One  way  to  conserve  energy  in  molecular  dynamics  runs  while  still  using  the  atom-by¬ 
atom  cutoff  method  is  to  set  the  atom-atom  potential  beyond  the  cutoff  distance  to  its  value  at 
the  cutoff  distance  (CUTOFF  II).  The  energy-force  pair  is  now  consistent  and  for  molecules 
like  water  where  the  forces  are  essentially  Coulombic  at  the  cutoff  distance  (and  the  total 
charge  on  each  molecule  is  zero)  the  energy  contribution  for  a  molecule-molecule  interaction 
conveniently  suns  to  zero  when  all  the  atom-atom  interactions  are  beyond  the  cutoff  distance. 
The  results  for  this  method  are  also  shown  in  Table  III  (CUTOFF  II).  Note  the  order  of  mag¬ 
nitude  reduction  in  the  standard  deviation  of  the  total  energy.  In  both  cutoff  methods,  waters 
have  a  tendency  to  "straddle”  the  cutoff  distance  boundary  in  such  a  way  as  to  reduce  repulsive 
and  increase  attractive  atom-atom  interactions.  For  CUTOFF  II,  this  has  a  much  smaller  effect 
as  the  potential  energy  for  any  atom-atom  interaction  changes  little  across  the  boundary.  For 
CUTOFF  I,  however,  each  atom-atom  potential  energy  function  is  truncated  to  zero  at  the 
cutoff  distance  which  causes  an  artificially  low  average  intermolecular  potential  energy. 

Another  method  which  might  seem  reasonable  to  try  in  order  to  create  a  consistent 
energy-force  system  for  molecular  dynamics  calculations  is  to  smooth  each  atom-atom  potential 
separately  to  zero  (AA  SMOOTH)  in  some  smoothing  region  and  then  take  the  derivative  to 
obtain  the  force.  Indeed  such  a  technique  might  be  useful  for  systems  where  the  value  of  the 
potential  at  the  cutoff  distance  is  near  zero.  For  water  this  not  the  case,  however,  and  AA 
SMOOTH  is  totally  useless  in  this  application.  For  our  test  run  we  smoothed  each  potential  to 
zero  from  0.85  nm  to  0.90  nm,  and  the  corresponding  force  was  calculated.  The  resulting 
energy  values  as  shown  in  Table  HI  differ  drastically  from  the  experimental  ones  due  mainly  to 
the  large  fluctuations  in  the  radial  distributions  near  the  cutoff  distance.  Large  forces  (20-30 
times  larger  than  for  the  unsmoothed  potential)  in  the  smoothing  region  cause  such  fluctuations 
and  are  a  result  of  the  steep  slope  of  the  potential  necessary  to  smooth  it  to  zero.  One  might 
view  this  effect  as  similar  to  smoothing  the  neglected  delta  function  force  term  of  the  CUTOFF 
I  system  over  0.05  nm.  Standard  deviations  are  not  given  for  AA  SMOOTH  because  the  total 
energy  of  the  system  continued  to  rise  over  the  course  of  the  run,  presumably  a  consequence  of 
the  large  forces  involved. 

The  technique  by  Andersen  and  Swope  (ANDERSEN)  which  smooths  each  entire  water- 
water  potential  to  zero  may  be  viewed  as  smoothing  the  delta  function  force  terms  of  a 
molecule-by-molecule  cutoff  or  equivalently  a  dipole-dipole  interaction  over  a  small  range.  0.05 
nm  in  our  case.  It  gives  the  best  energy  conservation  and  smallest  V,mra  and  Vymt  fluctuations 
as  shown  in  Table  III.  In  addition  its  waters  are  put  under  the  least  "stress”  as  measured  by 
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It  may  be  concluded  from  the  data  in  Table  III  that  there  are  several  advantages  to  using 
the  ANDERSEN  method.  It  should  also  be  noted  that  the  choice  of  method  of  handling  boun¬ 
dary  effects  significantly  influences  energy  calculations  with  differences  on  the  order  of 
10  kJ  mole-1. 

Table  III  also  contains  information  on  the  intramolecular  potential  energy,  the  bond 
lengths,  and  the  bond  angles  for  each  system.  V,m  may  be  seen  as  a  rough  index  to  the 
'stress*  each  water  molecule  is  experiencing.  For  an  ideal  gas  (IG),  i.e.  for  the  same 
intramolecular  potential  with  the  intermolecular  potential  turned  off,  the  intramolecular  poten¬ 
tial  energy  at  300  K  is  4.7  kJ  mole-1,  0.96  kJ  mole-1  above  the  (3/2 )k$T  value  of 
3.74  kJ  mole-1  one  would  expect  if  no  anharmonic  or  centrifugal  distortion  effects  were  found. 
In  the  liquid  state  each  oxygen -hydrogen  bond  on  the  average  is  stretched  and  each  HOH  angle 
on  the  average  is  reduced  below  the  equilibrium  value,  thus  increasing  V,m  above  the  ideal  gas 
level,  which  is  itself  above  the  harmonic  equipartition  value. 

For  a  non-rigid  calculation  at  295.4  K  using  the  LS  potential,  Rahman,  Stillinger,  and 
Lemberg72  report  4-  -  —34.8  kJ  mole-1,  while  we  calculate  for  WATTS  at  300  K 

*mm  +  Vim* "  —36.9  kJ  mole-1.  In  their  partitioning  between  Vim  and  Vlmir  they  assume, 
but  do  not  measure,  that  is  given  by  the  expected  undistorted  harmonic  oscillator  values, 
an  approximation  which  we  see  to  be  incorrect,  at  least  in  our  case,  due  to  anharmonicity,  cen¬ 
trifugal  distortion,  and  intermolecular  force  induced  molecular  distortion. 

The  experimental  value  to  which  the  calculated  intermolecular  potential  energies  should 
be  compared  deserves  some  discussion  as  two  significantly  different  numbers  are  quoted 
throughout  the  literature.  One  way  to  obtain  the  intermolecular  potential  energy  of  liquid  water 
is  to  equate  it  to  the  difference  in  energy  of  the  fluid  and  vapor  states.  This  is  calculated  by 
subtracting  PV  from  the  heat  of  vaporization  of  water  at  300  K.  Using  this  method,  Dashevsky 
and  Sarkisov104  obtain  for  the  intermolecular  potential  energy  from  experimental  data 
-41.0  kJ  mole-1  at  300  K,  and  -41.4  kJ  mole-1  at  298  K.  As  pointed  out  by  several 
workers,37*  l0°* 101* iaiM23  however,  the  bending  and  stretching  frequencies  of  mex  change 
upon  condensation,  and  this  difference  in  intramolecular  energy  must  be  acco&iMtu  for,  zn  *dl 
as  the  correction  for  conversion  of  free  to  hindered  translation  and  rotation.  Ranters  et  alK1 
estimate  a  correction  on  the  order  of  7.5  kJ  mole-1  which  would  lead  to  an  intermolecular 
potential  energy  of  -33.9  kJ  mole-1  for  298  K.  This  may  account  for  the  variation  in  experi¬ 
mental  intermolecular  potential  energy  quoted  in  the  literature,  as  some  workers  use  the 
corrected  -value  for  intermolecular  potential  energy  while  others  do  not.  It  should  also  be  men¬ 
tioned  that  the  definitions  used  of  Internal  energy”  and  "internal  potential  energy”  are  not 
always  well  spelled  out  or  consistent  among  different  authors,  making  comparisons  sometimes 
difficult  and  confusing. 

The  equivalent  experimental  value  for  total  energy  is  the  difference  in  energy 
between  liquid  water  at  300  K  and  ideal  noninteracting  water  vapor  at  0  K  with  no  zero  point 
vibrational  energy,  measuring  energies  from  the  bottom  of  the  potential  well  for  non-interacting 
molecules.  It  may  be  calculated  as  follows. 

£jook<«?)  -  £ok("W> 

m  lEj oox(fiff)  —  £ok(£*)1  +  [£o  *(£*)  —  £ox(vop)J  (4.13) 

—  —  Hok(*v)]  (4.14) 

because  E  H  for  liquid  water  and  ice.  Including  the  zero  point  vibrational  energy  gives  the 
results  shown  in  Table  IV. 


TABLE  IV.  Experimental  total  energy  (kJ  mole-1). 
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HoK.(ice)  —  Hot^wp)  -47.36* 

Vibrational  zero  point  energy  55.45* 
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*N.  Dorsey,  Properties  of  Ordinary  Water 
Substance,  (Hafner  Publishing  Co.,  New  York,  1968) 

*D.  Eisenberg  and  W.  Kauzman,  The  Structure  and 
Properties  of  Water,  (Oxford  University  Press,  New  York.  1969) 

Our  computed  EP"*  -  23.6  kJ  mole-1  and  experimentally  derived  E°*  -  21.5  kJ  mole-1 
total  energies  as  shown  in  Tables  II  and  IV  thus  agree  quite  well,  perhaps  better  than  expected 
in  light  of  the  possible  improvements  discussed  in  Section  V  below. 

F.  Heat  capacity 

The  energy  is  fixed  in  a  microcanonical  ensemble  while  the  temperature  as  computed 
from  Eq.  (2.6)  fluctuates  about  an  average  value.  Seven  distinct  water  configurations  with 
different  energies  are  created,  and  the  average  temperature  for  each  sample  is  calculated  over  at 
least  ten  picoseconds  of  running  time.  The  seven  points  are  plotted  on  the  energy-temperature 
graph  in  Fig.  7.  A  straight  line  is  fitted  to  the  points,  and  the  slope  is  calculated,  giving  the 
constant  volume  heat  capacity.  The  results  seen  in  Table  V  show  quite  good  agreement  with 
experiment  once  the  quantum  correction  is  added.  Note  that  the  calculated  value  would 
disagree  substantially  with  experiment  if  the  11.0  J  deg-Imole-1  intermolecular  quantum 
correction  for  hindered  rotational  and  translational  motion  had  been  omitted. 


TABLE  V.  Constant  volume  heat  capacity  (J  deg- 'mole-1). 


*D.  Eisenberg  ami  W.  Kauzman,  The  Structure  and 
Properties  of  Water,  (Oxford  University  Press,  New  York,  1969) 


V.  DISCUSSION  AND  CONCLUSION 

The  calculations  for  liquid  water  presented  here  are  designed  to  illustrate  the  quantum 
correction  of  classical  thermodynamic  quantities  and  not  to  provide  the  ultimate  in  accuracy  for 
those  thermodynamic  values.  Even  though  the  results  agree  well  with  experiment, 

-  23.6  vs  £**  -  21.5  kJ  mole-1,  and  C?"r  -  71.7  vs  CT  -  74.5  J  deg-1  mole-1,  it  is 
dear  that  these  classical  calculations  could  be  improved.  For  example,  it  can  be  argued  that  no 
potential  function  yet  exists  for  water  which  is  adequate  to  represent  both  the  inter  and 
intramolecular  motions  or  which  is  even  valid  in  an  effective  sense  for  ail  phases. 101  106  The 
WATTS  potential  function  which  we  use  in  this  example  calculation  is  no  exception,  having 
been  criticized 19  on  the  ground  that  radial  distribution  functions  calculated  from  it  do  not  agree 
with  experiment  It  is  unlikely,  as  we've  seen,  to  properly  account  for  the  change  in  vibrational 
frequendes10-37- 120,123  on  going  from  tire  gas  to  the  liquid  phase,  as  there  is  no  direct  coupling 
between  intermolecular  distances  and  the  intramolecular  part  of  the  potential.  The  reader  is 
referred  to  the  recent  paper  by  Reimers,  Watts,  and  Klein101  for  a  comparison  among  various 
existing  water  potentials  and  a  presentation  of  a  revised  Watts  potential.  The  potential  we  have 
used  is  dearly  only  an  effective35  71  molecule-molecule  potential,  as  it  omits  three  (ar.d  higher) 


molecule  effects35- 7*.  92*96. 124. 125  which  surely  must  exist.  In  addition,  one  could  make  a  more 
accurate  calculation  by  including  a  correction  6.16.37  for  the  tails  of  the  potential  beyond  the 
0.85  to  0.90  nm  region  at  which  we  feathered  the  potential  to  zero  or  one  could  try  other  long 
range  correction  techniques  such  as  Ewaid  sums.126  It  seems  dear  from  the  large  variations  in 
energy  among  different  choices  of  boundary  treatment  that  much  more  needs  to  be  learned 
about  the  effects  of  different  boundary  treatments  on  systems  with  long  range  potentials  and 
their  convergence  to  experimental  values.  Related  questions  have  been  raised  by  Pangali,  Rao, 
and  Berne65  with  respect  to  Monte  Carlo  calculations.  The  methodology  of  quantum  correction 
illustrated  here  Would  work  equally  well  with  any  or  all  of  the  improvements  mentioned  above 
to  the  classical  part  of  the  calculations. 

A  substantial  amount  of  calculation  is  needed  to  achieve  the  accuracy  illustrated  in  Fig.  7. 
The  long  simulation  time  to  achieve  a  stable  average  can  be  interpreted  in  terms  of  the  unusual 
"stickiness’  of  liquid  water.62*63  The  95  picoseconds  of  total  molecular  simulation  time  illus¬ 
trated  in  Fig.  7  required  190  hours  of  real  time  on  an  array  processor.107* 101  The  array  proces¬ 
sor  speed  is  approximately  35  times100  that  achieved  in  optimized  Fortran  on  a  DEC  VAX 
11/780  with  a  floating  point  accelerator,  and  judging  from  previously  reported  figures,62  5-10 
times  faster  than  a  rigid  water  calculation  on  an  IBM  360/91.  Our  2000  time  steps  per  hour 
when  scaled  for  number  of  particles  and  cut-off  radius  is  roughly  comparable  to  the  speed 
reported  by  Rapaport  and  Sdieraga*7- 110  for  their  array  processor  molecular  dynamics  calcula¬ 
tion  for  rigid  water,  taking  into  account  that  they  use  a  predictor-corrector  integrator,  while  we 
only  use  one  fence  evaluation  per  time  step.  The  total  number  of  atom-atom  force  evaluations 
is  1  x  10“  and  the  number  of  calculations  of  the  total  force  vector  on  an  atom  (summed  over 
all  its  pairwise  potential  interactions  with  other  atoms)  is  1  x  10*.  This  latter  figure  might  be 
roughly  compared  in  computational  effort  to  the  total  number  of  configurations  tried  in  a  simi¬ 
lar  Monte  Carlo  calculation,  i.e.  the  number  of  times  an  atom  is  moved  and  a  new  potential 
energy  is  calculated  as  a  sum  over  all  atomic  pairwise  interactions  with  other  atoms.  Since  each 
molecular  dynamics  atomic  force  evaluation  delivers  the  three  Cartesian  components  of  the 
atomic  force  vector  in  contrast  to  the  scalar  Monte  Carlo  computation  of  potential  energy,  it 
might  be  argued  that  the  proper  number  to  compare  to  an  equivalent  number  of  Monte  Carlo 
configuration  tries  is  3  x  10*.  One  might  also  compare  the  380  000  time  steps  to  an  equivalent 
number  of  Monte  Carlo  passes  through  all  variables.  It  has  been  argued  by  Rao,  Pangali,  and 
Berne62  that  one  Monte  Carlo  pan  for  rigid  water  can  be  compared  in  computational  effort  to 
one  molecular  dynamic  time  nap,  but  for  problems  accessible  to  Monte  Carlo  solution  that 
Monte  Carlo  may  be  several  times  more  efficient  in  terms  of  distance  moved  per  pass  versus 
per  molecular  dynamics  time  step. 

A  very  different  way  to  compute  dynamics  and  thermodynamic  quantities  which  may  in 
time  become  practical  would  be  a  quantum  force  classical  trajectory  approach127  in  which  at 
each  time  step  in  the  classical  trajectory  the  forces  (for  the  dynamics)  and  the  energy  (for  the 
thermodynamics)  are  computed  from  ab  Initio  quantum  mechanics. 

It  is  dear  from  these  results  that  one  can  and  should  take  into  account  quantum  correc¬ 
tions  in  testing  molecular  potential  energy  functions  against  experimental  thermodynamic  meas¬ 
urements.  In  particular,  the  intermoiecular  (hindered  translational  and  rotational)  motions  in 
strongly  assodated  liquids  can  lead  to  significant  errors  if  the  related  quantum  corrections  are 
neglected  in  thermodynamic  comparisons  with  experiment.  Consider,  for  example,  that  the 
intermoiecular  quantum  correction  to  energy  for  our  system  is  38  percent  of  the  kinetic  energy 
while  the  intermoiecular  quantum  correction  to  free  energy  is  20  percent  of  the  kinetic  energy. 
The  intermoiecular  quantum  correction  to  heat  capacity  is  15  percent  of  the  experimental  value 
while  the  intermoiecular  quantum  correction  for  entropy  is  10  percent  of  the  experimental 
value.109  Similarly,  motions  in  polymers  (which  can  themselves  be  affected  by  solvent  interac¬ 
tions)  may  also  need  thermodynamic  quantum  correction,  and  the  molecular  dynamics  method 
illustrated  here  also  can  be  applied  in  such  cases. 

An  interesting  aspect  of  this  quantum  correction  technique  is  that  after  the  dynamics 
(which  in  general  depend  upon  all  the  atoms)  are  computed,  the  quantum  corrections  may  be 
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calculated  atom  by  atom,  and  thus  the  quantum  effects  on  the  thermodynamic  variables  may  be 
considered  separately  for  different  elements,  different  chemical  environments  of  the  same  ele¬ 
ment,  different  types  of  molecules,  or  molecules  in  different  environments.  An  advantage  of 
the  basically  classical  molecular  dynamics  approach  to  thermodynamics  presented  here  is  the 
ability  to  visualize  and  understand  intuitively  the  classical  motions  and  frequencies  responsible 
for  thermodynamic  effects.  For  example,  one  can  understand  in  a  very  pictorial  way  the  domi¬ 
nance  of  the  water  quantum  corrections  by  the  hydrogen  atom  motions  as  illustrated  in  Figs.  3- 
6. 

This  technique  for  quantum  correcting  classical  thermodynamic  quantities  should  be  appli¬ 
cable  to  a  wide  variety  of  molecular  systems  including  polymers  such  as  proteins  and  nucleic 
adds,  liquids,  solutions  and  solids.  For  example  the  molecular  dynamics  method  could  be  used 
to  compute  and  quantum  correct  the  heat  capadty  of  biomolecules  in  solution,  a  quantity 
known  to  depend  on  molecular  conformation.  Thermodynamic  calculations  can  be  made 
involving  both  intermolecular  and  intermoiecular  degrees  of  freedom.  In  addition,  this 
approach  can  be  extended  to  treat  quasiequilibrium  cases,  such  as  the  calculation  of  thermo¬ 
dynamic  quantities  as  a  function  of  progress  along  a  chemical  reaction  coordinate  or  thermo¬ 
dynamic  quantities  for  molecules  in  special  surroundings  such  as  boundary  waters  near  a  pro¬ 
tein. 
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